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ABSTRACT 

Measured values of several variables are Incorporated 
into the calculation of snow water equivalent as 
"measured" from an aircraft by snow attenuation of 
terrestrial gamma radiation. Bayesian decision theory 
provides a "best" snow water equivalent measurement by 
taking into account the uncertainties in the individual 
measurement variables and "filtering” information about 
the measurement variables through prior notions of what 
the calculated variable (water equivalent) should be. 
Generalizations of principles involved in the applica- 
tion are presented in subsequent discussion. 

INTRODUCTION 

The importance of snow measurements for both water supply 
forecasting and for flood forecasting is well established. In 
the case of flood forecasting, the assessment of flood potential 
based on the water equivalent of snow on the ground provides a 
lead time for productive action which is not enjoyed during 
rainfall-caused floods. For example, during the record floods of 
the spring of 1969 in the upper Mississippi and middle Missouri 
River drainages, an estimated 137 million dollars in savings was 
realized through preventive action based on river forecasts. An 
estimated additional 97 million dollars savings was realized 
through the proper use of flood control structures based on fore- 
casts [Mondscheln, 1971]. Since snowwater equivalent is the 
primary data input to the forecast procedure, the importance of 
good snow measurements is readily seen. 

One very recent snow measurement method, the airborne gamma 
survey, appears to have great potential in plains areas such as 
the north-central United States. This method overcomes the dif- 
ficulties of nonrepresentative sample points and sparseness of 
samples now encountered in utilizing a ground-based point measure- 
ment network. Airborne measurements of gamma radiation coming 
from the soil in an area provide an excellent indication of the 
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amount of snow vater shielding the airborne gamma radiation 
detector from the radiation flux emanating from the ground. 

Either the flux in the entire energy range ("total count" 
method) , or the flux within a particular energy range ("spectral 
peak" method) can be used. The use of natural background radia- 
tion to measure snow water equivalent was first proposed in the 
Soviet Union [Kogan et al., 1965]. Dlmltriev et al. [1972] 
reported that operational airborne gamma snow surveys are now 
conducted over several million square kilometers of the U.S.S.R., 
mainly in regions where the ground hydrometeorological network 
is extremely sparse. Research in the method has been conducted 
in the United States since 1969 [Peck et al . , 1971; Burson and 
Fritzsche, 1972; Bissell, 1974; Larsen, 1975]. 

AIRBORNE GAMMA SNOW MEASUREMENT METHOD 

The total gamma radiation near the surface of the earth is 
dominated by flux of terrestrial origin. As seen in Figure 1, 
beyond an altitude of about 1500 meters radiation of cosmic 
origin begins to dominate. The spectral composition of terres- 
trial gamma radiation flux is shown in Figure 2. Prominent peaks 
in the spectrum indicate the presence of ^®K (1.46 MeV) , ^^^Bi 
(1.76 MeV, 2.20 MeV), and ^OSj.^ (2.62 MeV) in the lithosphere and, 
unfortunately, perhaps in the atmosphere as well. Radiation of 
atmospheric origin is unfortunate because it is a noise source in 
airborne gamma surveys. The isotope ^^^Bi is a subsequent decay 
product of the radioactive noble gas radon which, with a half-life 
of 3.6 days, can range far into the atmosphere from its origin 
uranium in the earth's crust. Consequently, much of the ^^^Bi 
decay radiation measured during snow surveys may not have been 
subjected to the hazard of attenuation by snow cover and is there- 
fore noise in the desired radiation signal. The isotope 208x1 is 
similarly the decay product of a noble radioactive gas (thoron), 
but is found in the atmosphere in only small quantities due to 
thoron's much shorter half-life (57 seconds). The primary snow 
measurement peaks are the spectral peak at 1.46 MeV and the 
208x1 spectral peak at 2.62 MeV. The radiation of atmospheric 
origin is of sufficient magnitude that the usefulness of the 
"total count" method is considerably reduced. The ^^K spectral 
peak also has an atmospheric contribution to due Compton-scatter- 
ing of 1.76 and 2.20 MeV photons do\m to the lower 1.46 MeV energy 
range, but a correction can be applied to "strip out" this effect 
by keying on the magnitude of the 1.76 MeV 214gx peak. Both the 
and 508xi peaks have cosmic components which can be "stripped 
out" by keying on the magnitude of the purely cosmic portion of 
the spectrum above 3 MeV. Another significant effect which can 
enter in is the variable attenuation of radiation within the soil 
Itself due to changing soil moisture, and a correction is also 
applied for this effect based on measured, estimated, or simulated 
soil moisture. 

The airborne "measured" value of water equivalent is given 
by the inverse of the relationship expressing the attenuation of 


410 



soil-originating spectral peak count rates as a function of water 
equivalent. An exponential form of this relation is given in 
equations (1) and (2): 
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effective water equivalent, including air 
blanket (g/cm'^); 

= attenuation coefficient of 1.46 MeV gamma 
rays in water (cm^/g); 

= attenuation coefficient of '^^°T1 2.62 MeV 
gamma rays in water (cm'^/g); 

= soil moisture (fraction of dry weight); 

® ratio of gamma attenuation cross section in 

water to cross section in air, taken to be unity; 
= expected value of count rate in decay in soil, 
with W=0 and S=0 (count s/ second) ! 

= expected value of count rate in 208^^ spectral 


peak due only to cos'll decay in soil with W=0 
and S=0 (count s/second ) ; 

The "stripping" equations (3) and (4) express the "pure" count 
rates Rg and Rjl as the actual spectral peak count rates 
obtained minus counts from noise sources: 
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where Cj^ = expected count rate in K spectral peak window; 

= expected value of a^craft background count rate 

contribution in spectral peak window; 

Cg = expected value of count rate in ^^^Bl spectral 
peak window; 

Bg = expected value of aircraft background count rate 
contribution in (1.76 MeV) spectral peak 

window; 

C’j’i = expected count rate in 208 j '2 spectral peak window; 

B^2 “ expected value of aircraft background count rate 

contribution in 208^]^ spectral peak window; 

Cq = expected value of count rate in cosmic (high energy) 
spectral window; 

^BK’^CK’^CT “ "stripping" coefficients, generally negative. 

Airborne water equivalent measurements are obtained by estimating 

Cjr, Cg, and Cq by the computed count rates Ng/t» 

Nj-j^/t, and Nj,/t, 

where Nj^ = counts in peak during time t; 

Ng = counts in 1.76 MeV ^^^Bi peak during time t; 
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208 

= counts In T1 peak during time t; 

Nj, = counts In cosmic window during time t. 

A linear combination of the and spectral peak measure- 

ments of a 15 cm. water equivalent will have an accuracy of some- 
where near one centimeter water equivalent [Bissell, 1974]. 

Operational use of the gamma snow measurements would require 
some quality control on the water equivalent calculation. Since 
automated calculations would be desirable, some mechanism which 
would automatically impose a constraint of "reasonableness" on 
the final calculated water equivalent would be in order. This is 
one area in which Bayesian decision theory seems to fit nicely. 

APPLICATION OF BAYESIAN DECISION THEORY 


Bayesian decision theory departs from classical statistics 
in that the distinction between parameters and random variables 
is not quite so sharp. It is recognized that initial parameter 
values may not actually be the true values, and some prior, or a 
priori, distribution for parameter values is allowed. 

Bayesian decision theory departs from classical estimation 
theory in that prior knowledge is allowed to play a role in the 
final decision. The difference between classical estimation 
theory and Bayesian decision theory is given in Figure 3. It is 
by introducing the prior distribution of the parameters that the 
"reasonableness" constraint is lnq>osed. A large variance on this 
distribution would allow the sample values a greater weight in the 
final estimate. A very small variance in the a priori distribu- 
tion would cause the estimation function to virtvially ignore an 
outrageous sample. 

In the case at hand, the random variables are the count 
values received in individual spectral windows. In vector form 
this is 


N = 


l”c J 


<5) 


The main parameter to be estimated is water equivalent, but other 
parameters must enter in as well because of their effect on KN. 
The expected values of the elements of N are given by 


= CRt 

= t{IU(W) + Bjj - BgjjCCg - Bg) - 

= t{R^*'^^ exp(-OjjW) + \ - exp(-oCgW) + Bg + G) 
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a 


B 


~ Bck^c^ (6) 

= rate of airborne radon contribution to ^^^Bi peak 
count; 

= expected value of count rate in ^^^Bi spectral peak 
due only to ^^^Bi decay in soil, with W=0, S=0; 

= attenuation coefficient of 1.76 MeV ^^^Bi gamma rays 
in water. 
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Looking at equations (6) through (9) , the items W (water equiva- 
lent) , S (soil moisture), G (representing the magnitude of ^Bi 
concentration in the atmosphere) , and C„ (cosmic flux) are the 
parameters. Thus the parameter vector is given by 


e = 



( 10 ) 


Now considering Q itself as a random variable, the desired result 
is to develop a posterior distribution for ® after the sample 
results are in [Llndgren, 1968, p. 234]. Let f (Q 1 g) be the con- 
ditional distribution of and g(6) be the a priori distribution 
of 0. The joint distribution of N and g is 


f(a.§) = f(n I 0)g(§) (11) 

The marginal (or absolute) distribution of N is given by 

hjjG) = Jf(a lg)g(£)dS (12) 

and finally the posterior density of S. conditioned on the observa- 
tion jj=a 

h(6| N) = f (a|S)g(8)/hj,(5) (13) 

If the loss function is taken to be (W-M) , then the Bayes risk 

(the posterior expected value of the loss function) is 

ffh[W-W(N)]^ = Bayes Risk (14) 

If we use the posterior expected value of W as the estimate of W, 

W = ^^(W) ■= |wh(al N)d§ (15) 


then this value minimizes the Bayes risk since the Bayes risk 
written above is a second moment, and the smallest second moment 
is that taken about the mean of a distribution [Llndgren, 1968, 
p. 239]. Hence is the Bayes estimate of the water equivalent. 

With the above background development, now for some simpli- 
fication. With a little (not unreasonable) bending, the problem 
can be fit into a multivariate normal framework, within which 
calculation of conditional means (l.e. the Bayes estimate - 
equation (5) - when a quadratic loss function - equation (4) is 
selected) is relatively simple. First it must be recognized that 
equations (6) through (9) express conditional means of N given 
Assume that g is multivariate normal. Analysis of the individual 
elements of g indicates that this is a reasonable assumption, 
with perhaps the normality of G falling the most in question. 
Also, W is an effective water equivalent which includes both the 
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air blanket and water blanket. Thus even for values of zero snow 
water equivalent , W Is not biunplng up against a physical boundary 
which would disallow negative values. 

Let the a priori (marginal) distribution of 0 be 

e^N(e*,lQ) (16) 


where 
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s 

G* 


and where 2(0 is known beforehand. Now (following Anderson, 1958, 
p. 28) let's put the random vector N and the parameter vector 0 
together Into a single augmented vector: 
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If X is multivariate normal. Its distribution Is designated as 


f (x)^iV()j,I) 


(18) 
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From (16), 



and 


(19) 


hi -h <“> 

When the distribution (18) is completely specified, calculation 
of the posterior expected value of 2 (our final goal) is a simple 
matter. First we must determine » £l 2 (which is the transpose 
of 521 ) £ 22 ' "^^sse are determined by looking at the 


414 



conditional mean and covariance of N based on 0. Continuing to 
follow the notation of Anderson, 


and 


( 21 ) 

- v(x^^^)][X^^^ - 

v*> ^ 

_1 def 

° ^22 ~ 52i5ii5i2 “ “22.1 


Let's see if we can identify some of the terms in (21) and (22) to 
be something known. First, recognize that equations (6), (7), (8) 
and (9) express conditional expectations. Assume that our prior 
knowledge of W and S is good enough that a fijst order expansion 
can be made about the prior expected values W and S . Then the 
exponential terms in eq\iations (6), (7), and (8) may be expanded. 
For example, 

R^^^exp(-Oj^W) ^ Rp^^exp(-ajjW*) K(S-S*) 

= * [ 1 + oqj(W-W ) 5- ] 

1-HCS 14KS 


(23) 

Thus by incorporating the linear expansions, adding and subtracting 
a few terms in the right places, and rearranging, we get 


^(5 10)4x1 = ^4x4 (§-«*)4xl+J'Jxl 


(24) 


A 

where the matrix A and vector § are constants. Comparing (24) 
and (21) we get the identification 


and 



i2i5n = ^ 


Therefore 

^21 “ 1 " 


(25) 

(26) 

(27) 


The elements of £22 i equation (22) are determined by physical 
reasoning. The matrix J 22 1 gives the covariance of the counts 
vector once g is specified! But once g is specified, the variabil- 
ity in the elements is due only to the random nature of nuclear 
counting- The nature of nuclear counting produces a Poisson random 
variable which, for sufficiently large expected values, is very 
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well approKimated by a normal distribution. This condition is 
satisfied in gamma snow surveys. Since the Poisson distribution 
has variance equal to the mean, we take the diagonal elements of 
^22 j to be the elements of the conditional mean vector of equation 
(21). Also, randomness of nuclear counting in one spectral window 
is independent of nuclear counting in another window during the 
same time. Thus the off-diagonal elements of ^ will be zero, 
and 


522.1 
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0 0 0 

V2(x<l)) 0 0 

0 V3(x(l)) 0 

0 0 


(28) 


Physical reasoning requires that the diagonal elements of E 22 1 
be positive. If the calculated values are not positive in any 
particular situation, a different method of computing water equiv- 
alent should be used. Finally then E 2 , is computed from equation 
(22) as 


-22 “ 522.1 52iSii5i2 


522.1 




(29) 


It should be noted that E 22 must be a positive definite matrix in 
order to be a marginal covariance matrix for N. Since J 22 1 
positive definite (diagonal, with positive elements), and since 
E is p.d. (by definition of marginal covariance matrix) it is 
sufficient to have matrix A with at least as many columns q 
as rows p, and be of rank p. This means that to compute J 22 
this fashion we must have at least as many parameter values (we 
have four) as count values (also 4), which condition is satisfied 
(Anderson, 1958, Appendix A, Theorem 1). The rank of the matrix 
A would have to be checked in each separate situation. 

Now to write the answer. The Bayes estimate of water equiv- 
alent under a quadratic loss function is the posterior expected 
value or, identically, the conditional mean based on the observed 
sample : 


A 

0 



= I X^^^) = + 2i2S22<5^^^- 

= ®* + 50 ^^^ 522.1 
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which Is seen simply to be a regression function on the count 
values. It is noted that not only is the posterior estimate of 
water equivalent given, but posterior estimates of soil moisture, 
atmospheric ^^^Bl concentrations, and mean cosmic flvix are given 
as well. Furthermore, the posterior covariance is given by a 
matrix computed by exchanging indices in equation (22). Thus 

the first Element in the first row of 2 is the posterior 
variance of the water equivalent measuremrat fj. 

THE KEY TO THE PROBLEM - ENCODING PRIOR KNOWLEDGE 

Looking back at the final product (equation (30)) of the 
previous section, it is seen that two major Inputs were necessary: 
(1) the matrix ^ which is a linearized version of the physical 
laws relating random variables and parameters, and (2) the vector 
0* and matrix Ig which quantify the description of prior 
knowledge. The whole key to the Bayesian approach is having 
enough prior knowledge to make all the gjnrations worthwhile. In 
the present application, the question is "where do we get our 
prior information?" The present formulation of the National 
Weather Service's Snow Accumulation and Ablation Model (Anderson, 
1973) will provide hydrological simulation of snow cover water 
equivalent W. Periodic ground measurements by cooperative 
observers where available assist in keeping the simulation "on 
track". It is anticipated that an associated soil moisture model 
will produce estimates of soil moisture S under the snow cover 
for input to the airborne gamma snow measurements. Development 
of the covariance matrix Eg first of all requires special calib- 
rations over a period of tme to quantify the performance of the 
hydrological simulations in updating previous measurements. The 
airborne 214gi concentrations certainly follow physical laws, but 
the inability to model these adequately means that in general some 
kind of a mean ^^^Bi level G* will have to be used in conjunction 
with a large "Ignorance factor" (entry on the diagonal of E^) for 
G. The cosmic flux is reasonably well defined as a function of 
season and barometric pressure. 

The most important point to be made in this paper is found 
in the above paragraph - the usefulness of simulation of physical 
processes to provide information which may be helpful in inter- 
preting signals received in the remote sensing process. This 
principle would seem to go beyond the airborne ganma snow survey 
application and extend to other remote observations of hydrological 
variables as well. The continuing effort to wring more and more 
information from remote sensing imagery would seem to require more 
and more sophistication as more and more automated interpretation 
processes are applied. A critical consideration is that progress 
must be made on a coordinated front so that remote sensing products 
are useful to the simulation practitioner and, similarly, that 
simulation results are available and useful to the remote sensing 
interpretation process. 
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Figure L. Reletlve contribution of terreetrlal 
and cosmic gamma flux with altlttide. 
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Figure 2. Spectral composition of terrestrial gamma radiation flux. 




r : 

I' : 


Figure 3. Classical estimation versus Bayesian decision 
theory to determine value of parameter 0. 
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